use "${clean}newspaper_analysis.dta", clear
keep if north == 1 | south == 1

// recoding Alexandria, VA to the South (technically part of DC until 1844)
replace south = 1 if city == "Alexandria"
replace north = 0 if city == "Alexandria"
replace state = "Virginia" if city == "Alexandria"

// dropping DC as a border area
drop if state == "District Of Columbia"

// standardizing variables
egen county_surban = std(county_urban_percent)
egen county_srugged = std(county_ruggedness)
egen county_spostoffice = std(county_postoffice_area)

// keeping the sample consistent
cap drop sample
gen sample = 1
foreach var of varlist $county_npcontrols {
	replace sample = . if `var' == .
}

// keeping only the north
keep if north == 1


// matrices to store the results for graphing
matrix define S = J(61, 3, 0)
*matrix define I = J(61, 3, 0)
*matrix define R = J(61, 3, 0)



// Loop
xtset pub_id
local row = 1

forvalues c = 1820 / 1880 {
	cap drop period
	cap drop time
	
	gen time = year - `c'
	gen period = time > 0
	
	qui xtreg singular c.time##i.period##ib1.lincoln_county ${county_npcontrols} if sample == 1, fe cluster(pub_id)
	
	*matrix R[`row',1] = `c'
	*matrix R[`row',2] = `e(r2_a)'
	*matrix R[`row',3] = `e(rmse)'	

	*qui lincom 1.period // intercept change at cutpoint
	*matrix I[`row',1] = `c'
	*matrix I[`row',2] = `r(estimate)'
	*matrix I[`row',3] = `r(p)'

	qui lincom 1.period#c.time // slope change at cutpoint
	matrix S[`row',1] = `c'
	matrix S[`row',2] = `r(estimate)'
	matrix S[`row',3] = `r(p)'
	
	local row = `row' + 1
}



// plotting the results
cap drop S1 S2 S3
svmat S
twoway ///
	(scatter S2 S1 if S3 < 0.05, msym(o) mcolor("${orange}")) ///
	(scatter S2 S1 if S3 > 0.05, msym(oh) mlcolor("${orange}")) ///
	, legend(ring(0) pos(1) order(1 "p < 0.05" 2 "n.s.") size(medsmall) row(1) region(lcolor(gs1))) ///
	xlab(1820 "1820" 1830 "1830" 1840 "1840" 1850 "1850" 1860 "1860" 1870 "1870" 1880 "1880") ///
	xtitle("Cutpoint Year", size(medsmall)) ///
	ytitle("Slope Change", size(medsmall)) ///
	subtitle("Lincoln Counties: Slope Changes at Each Year 1820 - 1880", box bexpand bcolor(gs2) color(white) size(medsmall)) 
graph export "${output}FigS6_np_lincoln_allslopes.pdf", as(pdf) replace









